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ABSTRACT 


Two  algorithms  are  presented  which  allow  for  the  unambiguous  resolution  of 
multiple  undersampled  frequency  components  in  a  signal.  Digital  signal  processing 
is  usually  governed  by  the  Nyquist  criterion  which  limits  the  amount  of  informa¬ 
tion  that  can  be  unambiguously  stored  and  recovered  digitally  to  a  spectral  width 
no  larger  than  half  the  sampling  frequency.  Both  algorithms  resolve  a  spectrum  be¬ 
yond  Nyquist  by  using  additional  information.  The  first  method  samples  a  signal 
more  than  once  using  a  different  sampling  frequency  each  time.  The  second  method 
utilizes  a  single  sampling  frequency  which  is  used  to  sample  both  the  signal  and  a 
band-limited  version  of  the  signal.  When  using  multiple  sampling  frequencies,  each 
sampling  frequency  yields  a  digital  sequence  which,  in  turn,  has  a  unique  spectrum 
when  the  Discrete  Fourier  Transform  (DFT)  is  applied.  The  bin  and  amplitude  infor¬ 
mation  from  each  of  the  resulting  undersampled  spectra  is  then  recombined  to  resolve 
the  original  spectrum.  In  like  manner,  when  using  a  single  sampling  frequency  the 
spectra  of  both  the  signal  and  its  band-limited  version  are  recombined  to  obtain  the 
solution.  Given  a  sampling  frequency,  both  algorithms  allow  for  the  unambiguous 
resolution  of  a  signal  with  a  spectral  width  at  least  twice  as  large  as  that  predicted 
by  Nyquist. 
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I.  INTRODUCTION 


A.  BACKGROUND 

For  a  periodic  discrete  signal  xp(nT )  with  a  period  of  N  samples  or  a  discrete 
signal  x(nT)  that  contains  exactly  N  samples,  a  transform  known  as  the  discrete 
Fourier  transform  (DFT)  is  defined.  The  DFT  transforms  the  N  time-domain  sam¬ 
ples  to  N  frequency-domain  coefficients  where  the  frequencies  are  2i rk/NT  where 
k.  =  0, 1, . . . ,  N  —  1  and  T  is  the  sampling  period.  The  DFT  coefficients  X(k)  are 
also  periodic  with  period  N.  For  actual  computations  using  a  digital  computer,  the 
continuous  signal  is  digitized,  the  number  of  samples  is  limited,  and  the  Fourier  rep¬ 
resentation  is  obtained  for  a  finite  number  of  frequency  values.  The  digitization  of  the 
signal  is  usually  governed  by  the  Nyquist  criterion  where  it  is  assumed  that  the  input 
signal  must  be  bandlimited  (0  <  /  <  fs/ 2)  before  going  into  the  analog-to-digital 
converter  (ADC).  The  Nyquist  theorem  however,  only  places  a  limitation  on  the  in¬ 
formation  that  can  be  derived  from  a  single  set  of  digitized  data  [3].  That  is,  a  single 
set  of  digitized  data  limits  subsequent  analysis  to  an  fs/2  bandwidth  unless  there  is 
additional  information  available.  With  additional  information,  the  frequency  compo¬ 
nents  /  >  fs/2,  which  appear  ambiguously  due  to  undersampling,  may  be  resolved 
provided  they  are  within  the  dynamic  range  of  the  algorithm  utilized. 

B.  UNDERSAMPLING 

There  are  several  advantages  to  an  undersampled  system  [4].  Among  these  are 
a  reduction  in  the  speed  requirements  on  the  digital  section  of  the  system,  a  relaxation 
on  the  analog  anti-aliasing  filter  requirements,  and  the  possibility  of  extending  the 
capabilities  of  existing  systems  with  relatively  minor  redesign.  Also,  a  good  deal 
of  power  and  money  can  be  saved  in  the  analog-to-digital  converter  section.  The 
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main  problem  in  identifying  the  frequencies  present  in  an  undersampled  signal  is  the 
resolution  of  the  ambiguities. 

There  has  been  a  great  deal  of  research  investigating  methods  to  either  de¬ 
tect  aliasing  or  to  recover  undersampled  signals.  Baier  and  Furst  [1]  have  proposed 
a  method  to  detect  aliased  frequency  components  by  utilizing  two  ADC’s  at  slightly 
different  sampling  rates.  The  frequency  components  occuring  in  the  resulting  spectra 
are  also  slighty  shifted  and  the  difference  in  the  two  spectra  allows  for  detection  of 
the  aliased  frequency  component.  Barbarossa  [2]  has  proposed  a  method  which  can 
recover  an  undersampled  signal  using  time-frequency  domain  analysis  and  computing 
the  Wigner-Ville  distribution  of  the  signal.  Rader  [10]  has  described  an  undersam¬ 
pling  technique  that  can  recover  periodic  signals  by  reconstructing  the  waveform  using 
a  set  of  trial  sampling  periods.  The  trial  period  which  yields  the  waveform  of  small¬ 
est  variation  is  then  considered  to  be  the  correct  period  and  the  resulting  waveform 
the  correct  waveform.  Zoltowski  and  Matthews  [13]  have  devised  an  algorithm  for 
unambiguous  estimation  of  both  frequency  and  phase  for  a  single  frequency  utilizing 
eigenvector  information.  To  come  more  in  line  with  real-time  wideband  processing, 
methods  based  on  the  use  of  phase  shift  information  to  resolve  the  ambiguities  in  a 
single  frequency  undersampled  signal  have  also  been  investigated  [12],  [11].  Finally, 
Martin  and  Rasmussen  [6]  have  proposed  a  testbed  for  evaluating  the  impact  that 
ADC  imperfections  would  have  towards  limiting  digital  implementations  of  under¬ 
sampling. 

C.  PRINCIPLE  CONTRIBUTION 

Current  research  in  the  area  of  undersampling  is  limited  in  two  principal  ways. 
First,  while  methods  exist  that  can  unambiguously  detect  a  single  frequency,  none 
of  the  research  proposes  a  dynamic  range  to  its  system  or  algorithm  for  detection. 
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Secondly,  no  algorithm  currently  exists  which  allows  for  the  unambiguous  detection 
of  a  signal  containing  a  spectrum  of  undersampled  frequencies. 

The  principal  contribution  of  this  thesis  is  to  resolve  both  problems.  This  thesis 
presents  algorithms  which  solve  for  an  undersampled  spectrum  of  frequencies.  The 
algorithms  are  based  upon  the  mathematical  properties  of  the  Symmetrical  Number 
System  (SNS)  [9]  and  the  relationship  between  digital  signals  and  the  SNS.  Each 
of  the  algorithms  developed  are  simple  linear  solutions  with  dynamic  ranges  which 
ensure  no  ambiguity  in  the  resolved  spectrum. 

D.  THESIS  OUTLINE 

This  thesis  presents  two  algorithms  which  unambiguously  resolve  an  under¬ 
sampled  spectrum  band-limited  beyond  Nyquist.  For  the  first  algorithm,  it  will  first 
be  shown  that  if  a  signal  is  sampled  at  two  different  sampling  frequencies  and  the  DFT 
is  applied  and  normalized,  two  unique  spectra  result.  Utilizing  the  bin  and  amplitude 
information  of  each  spectrum  as  well  as  knowledge  of  how  the  real  and  imaginary 
frequency  components  alias  into  DFT  bins,  it  is  shown  that  the  original  frequency 
components  may  be  unambiguously  resolved  utilizing  a  matrix  multiplication.  For 
the  second  algorithm,  it  will  be  shown  that  if  both  the  signal  and  a  band-limited 
version  of  the  signal  are  processed,  it  is  possible  to  resolve  a  spectrum  of  frequencies 
up  to  nearly  twice  the  Nyquist  rate.  The  band-limited  version  of  the  signal  allows  for 
the  unambiguous  resolution  of  frequencies  below  If  which  is  then  applied  to  solve  for 
frequencies  y  <  /  <  fs  —  1 

Since  the  properties  of  the  SNS  are  what  the  solution  methods  utimately  uti¬ 
lize,  Chapter  II  begins  with  a  detailed  look  at  the  Symmetrical  Number  System  (SNS) 
and  its  dynamic  range  for  unambiguously  detecting  a  tone  frequency  is  introduced 
and  proven.  Chapter  III  shows  that  the  DFT  encodes  frequency  information  in  a 
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format  that  is  exactly  the  same  as  the  SNS  which  means  that  a  digital  sequence  has 
the  same  properties  as  the  SNS.  Beginning  to  resolve  undersampled  frequencies,  the 
properties  of  the  SNS  allow  us  to  develop  the  architectures  of  Chapters  IV  and  V, 
in  which  two  and  three  channel  architectures  are  introduced  which  unambiguously 
resolve  a  single  frequency  within  the  SNS  dynamic  range.  Moving  towards  resolving 
multiple  frequencies,  Chapter  VI  shows  that  the  DFT  can  be  normalized  and  using 
the  properties  of  the  SNS,  a  method  is  demonstrated  which  resolves  two  undersam¬ 
pled  frequencies  in  a  signal.  With  the  SNS  and  normalization,  there  is  now  sufficient 
information  of  resolve  an  entire  spectrum  of  frequencies.  In  Chapter  VII  the  bin  loca¬ 
tions  and  normalized  amplitudes  are  used  to  set  up  linear  equations  which  solve  for  an 
undersampled  spectrum  with  frequencies  with  zero-phase  and  no  DFT  leakage.  Since 
real  world  signals  have  random  phase,  Chapter  VIII  removes  the  phase  constraint 
and  utlimately  presents  the  first  algorithm  which  resolves  frequencies  with  random 
phase  and  no  DFT  leakage  and  band- limited  at  twice  Nyquist.  Chapter  IX  presents 
the  second  algorithm  which  requires  only  one  sampling  frequency  but  has  the  same 
dynamice  range  as  the  first  algorithm.  Finally,  in  Chapter  X,  the  limitations  of  the 
two  algorithms  are  explored  and  a  possible  direction  for  resolving  the  limitation  is 
explored. 


4 


II.  THE  SYMMETRICAL  NUMBER  SYSTEM 


A.  DEFINING  THE  SNS  WAVEFORM 

In  the  SNS  preprocessing,  r  different  periodic  symmetrical  waveforms  are  used 
with  periods  based  on  r  different,  pairwise  relatively  prime,  integer  lengths  2, 
. . . ,  mr.  For  each  integer,  k,  an  r  dimensional  column  vector  Ak  is  formed  by  placing 
the  value  of  the  ith  waveform  at  k,  in  the  ith  position,  1  <  i  <  r.  It  is  then  desired  to 
find  out  the  largest  set  of  vectors  Aq, A\,...,Ak  that  are  distinct.  This  sequence  of 
k  +  1  vectors  forms  the  unambiguous  output  of  the  system.  This  number,  k  +  1,  is 
called  the  dynamic  range  of  the  system.  Dynamic  range  refers  to  the  fact  that  each 
vector  up  to  Ak  is  unique  in  the  same  way  that  frequencies  up  to  f8/ 2  are  unique  and 
unambiguous  when  sampled  and  conventional  methods  are  applied. 

The  definition  of  the  SNS  waveform  is  given  below. 

Definition  2.1  Let  m,  be  an  integer  greater  than  1.  For  an  integer  h  such  that 
0  <  h  <  m,  define 

ah  =  min{/i,m  —  h}  .  (1) 

This  function  is  extended  periodically  with  period  m.  That  is, 

®h+nm  =  a/i  (2) 

where  n  E  {0,  ±1,  ±2, . . .},  and  ah  is  called  a  symmetrical  residue  of  h  +  nm  modulo 
m. 


Let  am  be  the  row  vector  [ao,  01, . . . ,  am_  1].  For  m  odd 


0,1,..., 


m 


m, 


,...,2,1 


L  2  J  ’  L  2 

where  [x\  indicates  the  greatest  integer  less  than  or  equal  to  x.  For  m  even, 

mm,  ,  _  ' 

0, 1, -  1, ...  ,2, 1 


(3) 

(4) 
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Both  have  size  lxm  and  consist  of  the  symmetrical  residues  elements  ah,  0  <  h  <  m. 
This  shows  the  form  of  one  period  of  length  m.  From  this  definition  it  is  emphasized 
that  it  follows  for  any  integers,  h,  k,  ah  =  a,h+k  if  and  only  if  h  =  ±(h  +  k)( mod  m). 


B.  DYNAMIC  RANGE  OF  AN  SNS  SYSTEM 

Having  introduced  the  SNS  waveform,  it  is  useful  to  introduce  the  dynamic 
range  of  an  SNS  system  given  its  moduli. 


Theorem  2.1  Let  mi, . . . ,  mT  be  r  pairwise  relatively  prime  moduli,  and  let  Ao,  A\, 
A 2, ...  be  vectors  formed  by  the  symmetrical  number  system  given  in  Definition  2.1. 


a)  If  one  of  the  moduli  (m,\)  is  even,  then  the  dynamic  range  of  the  system  is 

M  =  min  j  II  +  II  7T?+  I  (5) 

[  1  e= 2  e=j+ i  J 

where  j  ranges  from  1  to  r  —  1  and  m,i2 ,  m,i3 , . . . ,  m,ir  range  over  all  permutations 
of  {2,3,...,  r}. 

b)  If  all  of  the  moduli  are  odd,  then  the  dynamic  range  of  the  system  is 

{l  j  i  r  I 

9  n  mie  +  9  II  mh  f  > 

1  £=i  z  e=j+ 1  J 

where  j  ranges  from  1  tor  —  1  and  ,  m,;2 , . . . ,  m,iT  range  over  all  permutations 
of  {1,2,...,  r}. 


1.  Example  of  the  SNS  Dynamic  Range 

Let  7779  =  4,  777,2  =  3,  and  7773  =  5.  We  must  minimize  the  set  of  values: 

f  777,1  777,!  777-1  ,  1  (,n  n  ,  ,  o7 

j—  +  777-2  '  777-3,  —777-2  +  777-3,  —  777-3  +  777, 2|  =  {17,  11,  Id}. 

The  dynamic  range  is  the  minimum  value  of  this  set,  11,  as  is  verified  in 
Table  2.1.  We  see  that  An  =  Alt  and  An  is  the  first  repetitive  vector.  In  other 
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Table  2.1:  Integer  Values  for  the  m\  —  4,  m-2  =  3,  and  m 3  =  5 

words,  Ao, . . . ,  A 10  is  a  set  of  11  distinct  vectors  and,  A0, . . . ,  An  is  not  a  set  of 
distinct  vectors.  This  shows  directly  that  the  dynamic  range  of  this  system  is  11.  A 
repeated  vector  is  called  an  ambiguity  of  the  system.  In  the  above  example  An  is  the 
first  ambiguity  of  the  system. 


C.  PROOF  OF  THE  SNS  DYNAMIC  RANGE 

Before  starting  the  proof  of  Theorem  2.1  it  is  necessary  to  express  the  general 
idea  that  will  be  used.  Suppose  that  are  r  positive  pairwise  relatively 

prime  integers.  Also  suppose  that 

’  4i  r  ai+k ' 

2  2 

.  ah  ah+k  . 

Ah=  .  .  =Ah+k  (7) 

.  ah  J  L  ah+k 


so  that  there  is  an  ambiguity  at  the  position  h+k,  (since  the  vector  is  indistinguishable 
from  that  at  position  h),  where  h  >  0  and  k  >  1. 

The  key  to  the  proof  is  that  fact  that  Ah  =  Ah+k  if  an  only  if  A,  =  ±(h  + 
k.)  (mod  m,j)  for  1  <  i  <  r.  Thus  the  proof  involves  systems  of  linear  congruences  and 
the  Chinese  Remainder  Theorem  [7].  The  goal  is  to  find  the  least  value  oi  h  +  k  that 
is  an  ambiguity.  In  general  all  permutations  of  subscripts  1, 2, . . . ,  r  must  be  worked, 
but  for  the  sake  of  a  clearer  presentation  the  subscripts  are  not  permuted. 
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1. 


Proof  of  Theorem  2.1a. 


Let  m  i  = 

2m  be  even. 

Case  I: 

h  =  (h  +  k)  (mod  m,), 

1  <i<j  (<  r) 

(8) 

h  =  —{h  +  k)  (mod  m,;), 

j  +  1  <  i  <  r 

(9) 

Case  II: 

h  =  —{h  +  k)  (mod  raj) 

1  <i<j 

(10) 

h  =  (h  +  k)  (mod  raj), 

j  +  1  <  i  <  r 

(11) 

Note  that  the  two  cases  are  opposite  in  terms  of  the  plus  or  minus  signs  in  the 
congruences.  These  would  not  be  separate  cases  except  for  the  fact  that  m,i  is  the 
even  modulus. 


a.  Case  I 

Suppose  congruences  (8)  and  (9)  are  true,  i.e.,  k  =  0  (mod  m* )  for 
1  <  i  <  j.  Since  the  moduli  are  relatively  prime  in  pairs,  k  =  0  (mod  Yil=\  mi)-  Note 
that  k  is  even  since  mi  is  even.  It  follows  that  k  =  aUUmi  for  some  a  G  {1,2,.. .} 
so 


k 


-  =  am  I] 


m7- 


(12) 


i= 2 


Continuing  on  for  the  moduli  m,j,  the  condition  i  >  j  + 1  is  examined  next.  From  (9) 


k 

h  = - (mod  m-i) 


(13) 


for  j  + 1  <  i  <  r.  Since  these  ml  are  pairwise  relatively  prime,  the  Chinese  remainder 
theorem  guarantees  that  there  is  a  unique  solution  h  (mod  Wi-j+\  mi)  to  this  system. 
Therefore,  there  will  be  exactly  Il£=imt  solutions  (mod  117=1  mi)-  Let  b'  be  the  least 
integer  such  that 

b'  m*  -  1  >  0  .  (14) 

i=j+ 1  z 
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Note  that  if  we  set 


r  k 

hi  =  b'  JJ  nti  -  -  , 
i=j+ 1  1 

(15) 

then 

T 

0  <  hi  <  JJ  m,i 
i—j+l 

(16) 

and 

k 

hi  =  --  (mod  mi) 

(17) 

for  j+1  <  i  <r.  This  hi  is  Resolution  (mod  n^+i m») 

.  All  solutions  (mod  nLi  md 

are  of  the  form 

T 

h  —  hi  +  b  JJ  m,i 

(18) 

i=j+ 1 


for  0  <  b  <  nLi  m,i.  Using  (15)  the  ambiguities,  h,  +  k,  will  have  the  form 

r  r  _  r  fc 

h  +  k  =  h\  +  b  JJ  m,i  +  k  =  b'  JJ  mi  +  b  JJ  m*  +  —  (19) 

i=7+l  i=j+l  *=i+l  ^ 

Using  (12) 

j  r 

h  +  k  =  am  JJ  m,i  +  b  m;  (20) 

i—2  i—j+l 

where  b  =  b'  +  b  €  {1, 2, . . .}  and  a  £  {1,2,...}. 

b.  Case  II 

Although  the  particulars  of  this  case  are  different  from  those  of  Case  I, 
the  line  of  argument  is  the  same,  and  the  result  is  exactly  the  same: 

j  T 

h  +  k  =  am  JJ  m,*  +  b  JJ  m,{  (21) 

i—2  i—j+l 

where  a  and  b  are  positive  integers. 

Now  the  two  cases  have  been  treated  and  it  is  evident  that  all  the  am¬ 
biguous  values  h  +  k,  have  the  form  (20)  or  (21)  and  these  forms  are  indistinguishable. 
The  least  number  of  this  form  (or  the  first  of  many  ambiguities)  is  found  by  setting 
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a  =  b  =  1.  Next  it  is  shown  that 


3  r 

JJ  m,i 

i= 2  i—j+ 1 


is  one  of  the  ambiguities. 

Suppose 


We  revert  to  Case  I  and  set 


m  n  m,  <  n  rrti 

i= 2  i=j+l 


a  =  1  in  (12)  so 


z  i=2 


Then,  in  (14)  b1  =  1,  so 


and 


*i=  n 


an,, - 


i—j+l 


hi  +  k  =  m  JJ  77?,j  +  JJ  7Tlj 
i=  2  i=i+l 

is  an  ambiguity.  That  is,  the  vector  corresponding  to  hi  is  the  same  as  the 
corresponding  to  hi  +  k. 

The  alternative  is  that 


r  3 

JJ  777/j  <  777  JJ  m,i  . 

i—j+\  i=2 

Now  we  use  Case  II  with  a  =  b  —  1.  We  can  set 


k 

2 


n 

i—j+l 


mi 


Then 


is  positive,  and 


5  A: 

hi  =  m  JJ  mi  -  - 

i=2  1 

3  r 

hi  +  k  =  mj[mi+  JJ  mi 

i— 2  i—j+l 


(22) 

(23) 

(24) 

(25) 

(26) 
vector 

(27) 

(28) 

(29) 

(30) 
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is  an  ambiguity. 


Recall  that 


’  a \  ' 

ah+k 

ah 

— 

ah+k 

l 

C3 

_ i 

.  ah+k  . 

(31) 


if  and  only  if  h  =  ±(h  +  k )  (mod  m*),  i  =  1 ,r.  There  are  2n  possibilities  for 
choice  of  plus  or  minus  signs  in  these  congruences.  For  j  fixed  we  have  looked  at  two 
possibilities  and  found  the  minimum  ambiguity.  All  others  can  be  found  by  varying 
j  and  by  permuting  the  order  of  the  moduli.  Thus 

M  —  min  j  m  ]1[  mit  1  (32) 

[  1=2  l=j+ 1  J 


is  the  minimum  ambiguity. 

2.  Proof  of  Theorem  2.1b. 

Let  TOi,...,mr  be  r  odd  pairwise  relatively  prime  natural  numbers.  Recall 


that 


'  ah  ' 

ah+k 

ah 

= 

ah+k 

.  “h  . 

.  ah+k  . 

(33) 


if  and  only  if  h  =  ±(h  +  k.)  (mod  m,;)  for  1  <  i  <  r.  Suppose  h  =  ( h  +  k )  (mod  m*)  for 


1  <  i  <  j  and  h  =  —{h  +  k)  (mod  for  j  +  1  <  i  <  r.  The  first  set  of  congruences 

implies  that  k  =  0  (mod  rrij),  1  <  i  <  j.  Thus 

j 

k  =  a  U  m,i  (34) 

2=1 

for  a  G  {1,  2, . . .}.  The  second  set  of  congruences  implies  that  2 h  =  —  k  (mod  m,;)  i.e., 


m 


-  1 


k,  (mod  m,i ) 


(35) 


for  j  +  1  <  i  <  r.  Again  the  proof  is  split  into  two  cases,  this  time  depending  on 
whether  k,  is  even  or  odd. 
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a.  Case  I 
Suppose  k  even.  Then 

h  =  ——  (mod  777, j)  (36) 

for  j  +  1  <  i  <  r.  Let  b'  be  the  least  integer  such  that 

hi  =  b'  rrii  -  ^  >  0  (37) 

i=j+ 1  1 

Note  b'  €  {1,2...}.  This  is  the  unique  solution  (mod  TlLj+i  mi)  t0  the  set  of 
congruences  (36).  Every  solution  to  the  original  problem  is  of  the  form 


h  =  b'  JJ  rrii  +  b  JJ 


mi~  2 


where 


0  <  b  <  rrii  . 


By  (34)  and  (38),  all  ambiguities  h,  +  k  that  arise  from  these  values  of  h,  are  of  the 


h  +  k  =  ^  |  m,i  (40) 

z  i= i  z  i=j+ 1 

where  b  =  2 (b'  +  b)  is  an  even  positive  integer.  Note  that  k  is  even  if  and  only  if  a  is 


b.  Case  II 

Suppose  that  k  is  odd.  As  in  Case  I,  all  ambiguities  h  +  k  that  arise 
from  these  values  of  h  are  of  the  form 

h  +  k,  =  ^  mi  +  \  II  mi  (41) 

z  i= i  1  i=j+ 1 


where  a  and  b  are  odd  positive  integers. 


Combining  Cases  I  and  II,  note  that  the  smallest  value  of  (40),  (41)  is 
found  when  a  =  b  =  1.  Furthermore,  when  a  =  b  =  1,  and  n(=i  mi  <  Ui=j+imi  this 


12 


is  an  ambiguity:  Let 


h  =  \  f[  rni-^Urm  (42) 

z  i=j+ 1  z  i= 1 

and  /c  =  ]li=imi-  It  is  easily  checked  that  h  =  h  +  k(mod  m*)  for  1  <  i  <  j  and 
/i  =  —  (/?,+  fc)  (mod  m,i)  for  j  +  1  <  i  <  r.  When  n[=j+i  mi  <  lli=i  mi  we  must  solve 
the  permuted  problem:  h  =  /i+A:(mod  mi)  for  j  +  l  <  i  <  r  and  h,  =  —  (A+A:)(mod  m,) 
for  1  <  i  <  j.  This  time  k  =  IILj+i^i,  and  h  =  (l/2)rii=imi  _  (1/2)  IIi=j+i  mi 
give  the  ambiguity  at  h  +  k.  Since  the  solution  to  the  problem  involves  looking  at 
all  permutations  this  gives  the  minimum  ambiguity.  This  completes  the  proof  of 
Theorem  2.1b. 

Theorem  2.1  provides  the  basis  for  predicting  the  dynamic  range  of 
any  system  which  encodes  information  in  a  format  identical  to  the  SNS.  In  fact,  in 
the  next  chapter  it  is  shown  that  the  process  of  digital  sampling  encodes  frequency 
information  into  the  SNS  format  which  allows  the  DFT  to  be  treated  identically  to 
the  SNS. 
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III.  DFT  AND  SNS  RELATIONSHIP 


A.  DIGITAL  SAMPLING  AND  ALIASING 

Consider  a  single  frequency  signal  sampled  at  two  different  sampling  frequen¬ 
cies.  Digital  uniform  sampling  of  an  analog  waveform  with  frequency  /  produces  a 
discrete  waveform  which  is  symmetrical  about  the  sampling  frequency,  fsJ 2.  Assume 
for  this  system  that  the  two  sampling  frequencies  are  fs i  =  10  and  fs 2  =  11.  After 
sampling,  an  analog  input  signal,  x(t),  becomes  a  discrete  sequence,  x{nT).  This 
periodic  sequence  has  a  digital  frequency  given  by  u>  =  2n(f  / fs).  A  signal  with  dig¬ 
ital  frequency  0  <  u>  <  -k  is  indistinguishable  from  a  signal  with  digital  frequency 
nn  <uj  <  (n  +  l)7r,  n  =  1, 2,  3, . . .,  an  effect  known  as  aliasing. 

B.  DIGITAL  FREQUENCY  MAPPING 

The  digital  frequency  of  a  sampled  sinusoid  can  be  mapped  into  the  z-domain 
as  shown  in  Figure  3.1a.  For  simplicity  assume  a  sinusoid  x(t )  =  2  cos  2a  ft  and  after 
sampling 

x{n)  =  2  cos  urn  =  ejijn  +  e~jujn  .  (43) 

If  /  =  fs/ 4,  this  corresponds  to  u  =  7r/2.  If  /  =  fs/ 2,  w  =  tt.  Since  the  signal  is 
real,  the  signal  appears  in  complex  conjugate  pairs  on  the  z-plane.  For  frequencies 
between  fs/2  and  /s,  the  frequencies  map  back  to  their  conjugate  on  the  upper  half 
of  the  complex  plane.  If  the  frequency  is  increased  beyond  fs ,  a  full  trip  is  made 
around  the  unit  circle  and  the  mapping  repeats.  Figure  3.1b  illustrates  the  mapping 
with  each  triangle  representing  a  full  rotation  around  the  unit  circle  in  the  z-plane. 
The  abscissa  represents  the  input  analog  frequency  while  the  ordinate  represents  the 
digital  frequency  mapping.  Note  that  an  infinite  number  of  analog  frequencies  will 
map  into  each  digital  frequency  0  <  to  <  tt. 


15 


Imag 


Figure  3.1:  a)Z-plane  mapping  of  an  input  analog  signal,  b)  sampled 
frequency  output. 


C.  DFT  BINS  AND  THE  SNS 


Recall  that  the  DFT  is  given  by: 

N- 1 

X(k)=J2  x(n)e~j^nk/N)  k  =  0, ...  ,N  —  1  .  (44) 

n=0 


Application  of  the  DFT  to  x(n )  yields  a  discrete  spectrum  where  X(k)  is  the  energy 
contained  in  the  signal  at  each  digital  frequency  u  =  2nk/N .  The  discrete  spectrum 
X  ( k )  has  N  indices  with  the  digital  frequency  of  each  index  given  by: 


0,  2tt 


N  ’ 


,2tt 


(N/2)  n  (N/ 2  +  1) 


N 


27 r 


TV 


27T 


(iv  —  2)  27r_(iv-i) 


A 


N 


for  AT  even 


(45) 
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and 


1  (TV 

°-  — jf 


l)/2  2JN  + 1)/2 


TV 


TV 


TV 


for  TV  odd . 
(46) 


The  analog  frequency  corresponding  to  each  index  is  obtained  by  multiplying  each 
value  by  fs.  Since  signals  with  digital  frequencies  in  the  range  tt  <  lu  <  2tt  are  indis¬ 
tinguishable  from  signals  with  digital  frequencies  0  <  u>  <  tt,  the  digital  frequency  of 
each  index  can  also  be  written  as: 


0,  2t r 


TV’ 


TV 


TV 


■  2wb  2*i7 


for  TV  even  ,  (47) 


and 


0,  2t r 


TV’ 


2k 


[N/2]  ITV/2J 


TV 


2ft- 


N 


_  2  1 
’  2nN’  2*N 


for  N  odd 


(48) 


where  |_xj  represents  the  greatest  integer  less  than  or  equal  to  x.  More  simply,  the 
spectrum  X(k)  resolves  into  N  integer  indices  and  incoming  signals  will  map  into 
unique  bins: 


AT  TV 
2  ’  2 


for  TV  even  , 


(49) 


and 


,  TV,  ,  TV  ,  ' 


for  TV  odd  . 


(50) 


That  is,  since  frequency  indices  greater  than  TV/2  are  redundant  for  real  signals, 
the  highest  unaliased  frequency  that  can  be  observed  corresponds  to  the  TV/2  index. 
From  the  above  discussion,  it  is  clear  that  the  DFT  maps  real  signals  naturally  into 
the  symmetrical  number  system  [9].  In  this  case,  the  modulus  described  in  [9]  and 
in  Chapter  II  is  the  sampling  frequency,  fs.  The  number  of  indices  TV  is  given  by 
TV  =  / sTl ,  where  Tl  is  the  total  sampling  time. 


D.  ILLUSTRATING  THE  RELATIONSHIP 

Figure  3.2  illustrates  the  DFT  mapping  for  two  channels  where  fs i  =  10  and 
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fS2  =  11  for  input  frequencies  /  =  0  to  23.  In  this  case,  Tl  =  1  so  N\  =  10  and 
iV2  =  11.  In  the  figure,  the  abscissa  corresponds  to  the  incoming  frequency  while 
the  ordinate  corresponds  to  the  bin  the  signal  is  resolved  into.  Table  3.1  displays  the 
input  frequency  and  the  resulting  DFT  bin  for  each  sampling  frequency.  Note  that 
frequencies  resolve  as  described  in  (49)  and  (50).  By  considering  both  channels,  it  is 
possible  to  unambiguously  resolve  signal  frequencies  in  the  dynamic  range  determined 
by  the  SNS,  (0  <  /  <  15). 

Simply  knowing  that  the  DFT  is  in  the  SNS  format  and  what  the  dynamic 
range  is  for  any  given  sampling  frequencies  would  be  useless  unless  the  information 
can  be  used  to  resolve  undersampled  frequencies  within  the  dynamic  range.  The  next 
two  chapters  detail  two  and  three  channel  architectures  as  well  as  solution  methods 
which  will  resolve  a  single  undersampled  frequency  band-limited  within  the  dynamic 
range  of  the  SNS  system. 
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Input 

Frequency 

DFT  Bins 

7 

fs  =  10 

fs  =  11 

0 

0 

0 

i 

1 

1 

2 

2 

2 

3 

3 

3 

4 

4 

4 

5 

5 

5 

6 

4 

5 

7 

3 

4 

8 

2 

3 

9 

1 

2 

10 

0 

1 

11 

1 

0 

12 

2 

1 

13 

3 

2 

14 

4 

3 

15 

5 

4 

16 

4 

5 

17 

3 

5 

18 

2 

4 

19 

1 

3 

20 

0 

2 

21 

1 

1 

22 

2 

0 

23 

3 

1 

Table  3.1:  Input  Frequency  and  Resulting  DFT  Bins  for  2  Channel  Exam 
pie 
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IV.  RESOLVING  THE  TWO  CHANNEL  CASE 


A.  THE  GENERAL  SETUP 

Figure  4.1  shows  the  block  diagram  of  a  two-channel  receiver  architecture  to 
determine  a  single  frequency  /.  In  this  architecture  the  ADC  sampling  frequencies 
fs i  and  fS2  are  relatively  prime.  The  DFT  outputs  are  thresholded  to  detect  the 
frequency  bins.  The  frequency  bins  ai  and  a2  are  then  used  by  the  SNS-to-decimal 
algorithm  to  determine  the  frequency  of  the  input  signal.  Let  mi  =  fs  1  and  m2  =  /s2 
and  suppose  that  the  incoming  frequency,  /  (unknown)  lies  within  the  dynamic  range 
M  of  the  SNS  system  (5),  (6).  Thus,  /  =  ±ai(mod  mi)  and  /  =  ±a2(mod  m2).  For 
each  of  these  congruences  either  the  plus  or  the  minus  is  correct,  but  it  is  not  known 
which.  Thus,  there  are  four  sets  of  two  equations: 


i) 

/ 

= 

ai(mod  mi), 

/ 

= 

a2(mod  m2), 

ii) 

/ 

= 

ai(mod  m-i ), 

/ 

-a2(mod  m2), 

iii) 

/ 

= 

— ai(mod  mi), 

/ 

= 

-a2(mod  m2), 

iv) 

/ 

= 

— ai(mod  m-i), 

/ 

= 

a2(mod  m2)  . 

The  Chinese  remainder  theorem  [7]  guarantees  that  each  of  these  has  a  unique  solution 
modulo  mim2,  and  Theorem  2.1  guarantees  that  exactly  one  of  these  solutions  lies 
within  the  dynamic  range  of  the  system  and  this  is  the  value  of  /.  In  fact,  it  is  only 
necessary  to  solve  i)  and  ii),  at  most,  because  the  solutions  to  iii)  and  iv)  are  the 
negatives  of  the  solutions  to  i)  and  ii),  respectively. 

Recall  that  in  the  standard  statement  of  the  Chinese  remainder  theorem  we 
wish  to  solve  for  /,  where  /  =  Oj(mod  m,),  1  <  i  <  r,  and  the  m;  are  pairwise 
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Figure  4.1:  Block  diagram  of  a  two  channel  receiver  architecture  to  deter¬ 
mine  a  single  frequency  / 

relatively  prime.  The  theorem  states  that  there  is  a  unique  solution  modulo  M  — 
7rcim2  •  •  ■  m,r.  A  standard  method  of  solution  is  to  find  integers  b{  such  that 

=  l(mod  m,i)  ,  (51) 

777-j 

i  =  1, 2,  •  -  • ,  7 ■,  in  which  case  the  solution 

f  =  Mb\a  i/m,i  +  Mb2(i2/'m'2  ■■■  T  M  brar/Tn,r  (mod  M )  .  (52) 

Returning  to  the  two  channel  case,  note  that  the  values  of  b\  and  62  depend 
only  on  777-1  an<f  and  not  at  all  on  ±ai  or  ±02-  Thus  it  may  be  assumed  that 
the  constants  cx  =  m,2bi(=  Mb\/m{)  and  c2  =  777-162  are  known,  and  that  the  SNS- 
to-decimal  algorithm  only  needs  to  evaluate  ±ciai  ±  c2a2  modulo  M(=  777-1777-2),  and 
pick  the  one  value  that  lies  within  the  dynamic  range. 
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B.  SET-UP  FOR  A  PARTICULAR  TWO  CHANNEL  RECEIVER 


Let  m i  =  330  and  m2  =  337.  Then 

bi  =  -47  solves  3376x  =  1  (mod  330)  ,  (53) 

and 

b2  =  48  solves  330^2  =  1  (mod  337)  .  (54) 

Thus,  ci  —  m2bi  =  -15,839,  and  C2  =  m,ib2  =  15,840.  Also,  M  —  111,210,  and,  by 
Theorem  2.1,  the  dynamic  range  of  the  system  is  m,i/2+m,2  =  502.  These  coefficients 
are  hardwired  and  do  not  need  to  be  recomputed. 

C.  SOLVING  A  PARTICULAR  CASE  -  THE  WORK  OF  THE  ALGO¬ 
RITHM 

Suppose  that  the  frequency  of  the  incoming  signal  is  /,  0  <  /  <  501,  and 
ai  =  59  and  a2  =  66.  By  Theorem  2.1,  /  is  the  unique  solution  to  ±(—15, 839)(59)  ± 
(15,840)(66)  modulo  111,210  that  lies  in  the  interval  [0,501].  First  try  (  —  15,839)  • 
59  +  15,840  •  66  =  110,939;  this  is  congruent  to  —271  modulo  111,210,  so  case  iii) 
solves  /  =  271. 

It  is  emphasized  that  this  is  the  only  computation  the  algorithm  needs  to 
compute. 

D.  SPECIAL  CASES 

In  special  cases  there  are  even  quicker  solutions  than  given  in  section  C.  For 
example,  suppose  m\  =  2 p  and  m2  =  2p  +  1.  By  Theorem  2.1  the  dynamic  range  of 
the  system  is  m\j 2  +  m,2  =  3p  +  1.  Frequencies  within  the  dynamic  range  will  fall 
into  bins  as  follows: 

(55) 


p+1  terms 


p  terms 
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p  terms 


The  top  bin  represents  the  sampling  frequency  of  m\  =  2 p  and  the  bottom  m2  =  2p+l. 
Now  suppose  the  incoming  frequency,  /,  is  resolved  into  bins  ai  and  02,  respectively. 
From  the  above  it  is  clear  that: 

{CL\  if  Ul  =  &2 

mi  —  01  if  ai  =  a2  —  1  (56) 

m,i  +  ai  if  01  =  02  +  1 

This  makes  the  solution  much  more  straightforward. 
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Figure  5.1:  Block  diagram  of  a  three  channel  receiver  architecture  to  de¬ 
termine  a  single  frequency  / 

or 

/  =  ±252  ±  70  ±  120  (mod  210)  ,  (58) 

but  252  (mod  210)  =  42  (mod  210)  so: 

/  =  ±42  ±  70  ±  120  (mod  210)  .  (59) 


First, 

/  =  42  +  70  +  120  =  232  =  22(mod  210)  ,  (60) 

a  value  that  can  be  discarded,  because  it  and  its  negative  are  out  of  the  dynamic 
range. 

Second, 

/  =  42  +  70  -  120  =  -8(mod  210)  .  (61) 

Although  -8  is  out  of  range,  the  negative,  /  =  8,  is  in  the  dynamic  range,  so  that 
/  =  8  is  the  correct  frequency  value. 


26 


VI.  RESOLUTION  OF  TWO  UNDERSAMPLED 

FREQUENCIES 


A.  DFT  AMPLITUDE  NORMALIZATION 

Consider  a  unit  amplitude  signal  containing  a  single  integer  frequency.  After 
sampling  and  forming  an  N  point  DFT  X(k),  the  kth  bin  corresponding  to  the  digital 
frequency  of  the  signal  will  have  an  amplitude  A(k)  given  by: 


A(k)  = 


N  if  k  =  0,  or  Neven/2 

N/2  elsewhere 


(62) 


Dividing  X(k.)  by  A(k.)  will  produce  a  normalized  spectrum. 

Figure  6.1  shows  two  DFTs  of  an  undersampled  signal  containing  two  fre¬ 
quency  components  with  unknown  amplitude  and  frequency.  Figure  6.1a  displays  the 
frequency  spectrum  for  /sl  =  10  while  Figure  6.1b  corresponds  to  fS2  =  11.  The  two 
frequencies  are  known  to  lie  within  the  dynamic  range  of  the  SNS  system  which  in 
this  case  corresponds  to  /  <  15.  It  is  difficult  to  determine  which  two  SNS  pairs  are 
valid  by  looking  at  the  spectra  of  Figure  6.1.  That  is,  either  the  bin  pairs  \  ,  [|1  are 
correct  or  [3] ,  [2]  are  correct. 


B.  BENEFITS  OF  NORMALIZATION 

The  benefits  of  normalizing  the  DFT  are  shown  in  Figures  6.2a  and  6.2b  where 
each  of  the  original  spectra  have  been  divided  by  A\{k)  and  A2{k)  respectively.  In 
this  case  Ni  =  10  and  N2  =  11.  It  is  now  easy  to  see  which  magnitudes  “match-up” 
and  that  \  ,  3  is  the  correct  set.  Using  the  special  case  of  two  SNS  moduli  presented 
in  Chapter  IV,  mi  =  2 p  and  m2  =  2p  +  1  where  p  =  5,  it  is  determined  that  [j  ,  \ 
corresponds  to  f=9  and  f=14  with  amplitudes  0.7514  and  0.6515  respectively. 


Theorem  6.1  Let  fs\  and  fs2  be  any  two  relatively  prime  sampling  frequencies.  Con¬ 
sidering  a  signal  containing  two  unknown  integer  frequencies  with  distinct  amplitudes, 
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the  N  -point  normalized  DFTs  reveal  two  unique  bin  pairs  which,  in  turn,  reveal  the 
two  frequencies  provided  they  are  within  the  dynamic  range  of  the  SNS  system. 


C.  SPECIAL  CASE 

The  solution  has  a  special  case  when  the  two  frequencies  alias  to  the  same  bin 
for  one  of  the  sampling  frequencies.  For  example,  incoming  frequencies  /  =  3  and 
f  —  13  would  both  alias  into  bin  3  for  fs i  =  10  but  into  bins  3  and  2  for  fs 2  =  11. 
Normalization  of  the  DFT  will  still  produce  accurate  results  as  the  magnitudes  at 
bins  2  and  3  for  fs 2  =  11  will  add  up  to  the  magnitude  of  bin  2  for  fs±  —  10  and  the 
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b) 


Figure  6.2:  Normalized  DFT  for  two  unknown  frequencies  for  a)  fs i  =  10 
and  b)  /s2  =  11. 

bin  pairs  3],^]  will  yield  the  correct  frequencies.  When  two  frequencies  are  received 
with  identical  amplitudes,  Theorem  6.1  does  not  apply. 

Having  shown  that  normalization  allows  us  to  “see”  which  bin  pairs  “match” , 
normalization  can  also  be  used  to  measure  how  much  energy  went  into  each  bin.  If 
two  frequencies  land  in  the  same  bin,  the  normalized  energy  should  essentially  add 
together  into  that  bin.  This  information  along  with  bin  location  is  used  in  the  next 
chapter  to  solve  for  an  entire  spectrum  of  frequencies. 
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VII.  INTRODUCTION  TO  RESOLVING  MULTIPLE 

FREQUENCIES 


A.  SETTING  UP  THE  SOLUTION 


The  normalization  of  the  DFT  amplitude  introduced  in  the  previous  chapter 
leads  to  the  consideration  of  a  method  to  resolve  any  number  of  zero-phase  frequencies 
which  have  no  DFT  leakage.  By  zero-phase  frequencies  with  no  leakage  it  is  meant 
that  the  DFT  values  will  be  all  real  and  that  each  resolves  its  entire  energy  into  one 
DFT  bin  with  no  leakage. 

Consider  a  bandlimited  signal  with  frequency  components  that  can  range  from 
0  to  10  Hz.  The  amplitude  and  frequency  of  any  component  is  unknown  except  that 
they  are  integer  frequencies.  Let  x0,  aq,  •  •  • ,  x10  represent  the  amplitudes  of  each 
frequency  component  from  f=0  Hz  to  f=10  Hz.  Assume  that  the  signal  is  sampled 
at  /si =10  Hz  and  fs 2=11  Hz  and  that  the  DFT  is  applied  and  normalized.  Two 
normalized  spectra  will  be  formed,  each  containing  undersampling  ambiguities  for 
/  >  5.  Let  the  notation  b^m  represent  the  value  of  the  normalized  DFT  coefficient 
for  bin  i,  fs  =  m.  From  knowledge  of  aliasing  it  is  clear  that: 


In  other  words, 


£>i,io 

= 

X\  +  £9 

£>2,10 

= 

X2  +  X8 

£>3,10 

= 

Xz  +  x7 

£>4,10 

= 

£4  +  Xq 

£>5,10 

= 

X5 

£>5,11 

= 

X5  +  XQ 

£>4,11 

=r 

X4  +  £7 

£>3,11 

= 

xz  +  x8 

£>2,11 

— 

X2  +  Xq 

£>i,u 

— 

Xl  +  £10 

£>o,n 

= 

£0 

for  fsi  =  10,  /  =  1  and  /  =  9  alias  into  bin  1,  while 


(63) 


/  =  2 


and  /  =  8  alias  into  bin  2,  and  so  on.  This  forms  a  set  of  linear  equations  that  can 


31 


be  written  as: 


'  0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

0  ' 

x0 

1 

0- 

T-* 

0 

_ 1 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

®1 

62,10 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

x2 

h,io 

0 

0 

0 

0 

1 

0 

1 

0 

0 

0 

0 

x3 

&4,10 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

X4 

65,10 

0 

0 

0 

0 

0 

1 

1 

0 

0 

0 

0 

x5 

= 

65,11 

0 

0 

0 

0 

1 

0 

0 

1 

0 

0 

0 

x6 

fr4.ll 

0 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

x7 

fr3.11 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

x8 

fr2.11 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

Xg 

fri.n 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0  . 

.  ®10  . 

.  fro, 11  . 

(64) 


This  is  of  the  form  Ax  =  b  where  x  represents  the  actual  amplitudes  of  each  frequency 
component  and  b  represents  the  total  amplitude  that  is  aliased  into  each  bin  for  the 
two  sampling  frequencies.  From  (63)  it  is  easily  seen  that  the  values  x0  through  xio 
may  be  uniquely  solved  for.  To  begin  with  xo  =  6o,n  and  X5  =  65,10-  Then  the  values 
xq,  X4,  xj,  etc.,  may  be  solved  for  successively.  Thus,  A  is  nonsingular  and  x  =  A~lb: 


x0 

'  0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 ' 

frl,10 

Xi 

1 

1 

1 

1 

1 

-1 

-1 

-1 

-1 

0 

0 

fr2,10 

0 

1 

1 

1 

1 

-1 

-1 

-1 

0 

0 

0 

*>3,10 

0 

0 

1 

1 

1 

-1 

-1 

0 

0 

0 

0 

fr4,10 

£4 

0 

0 

0 

1 

1 

-1 

0 

0 

0 

0 

0 

fr5,10 

*5 

= 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

fr5,ll 

Xq 

0 

0 

0 

0 

-1 

1 

0 

0 

0 

0 

0 

fr4,ll 

£7 

0 

0 

0 

-1 

-1 

1 

1 

0 

0 

0 

0 

fr3,U 

*8 

0 

0 

-1 

-1 

-1 

1 

1 

1 

0 

0 

0 

fr2,ll 

£9 

0 

-1 

-1 

-1 

-1 

1 

1 

1 

1 

0 

0 

frl.ll 

.  *10  . 

_  -1 

-1 

-1 

-1 

-1 

1 

1 

1 

1 

1 

0 . 

.  fro, 11  . 

(65) 


From  (65),  it  is  now  possible  to  reconstruct  the  original  spectrum  based  upon 
the  normalized  bin  values  for  the  two  undersampled  spectra.  Where  before  only 
frequencies  from  0  to  5  Hz  could  be  recovered,  65  allows  integer  frequencies  from  0  to 
10  Hz  to  be  unambiguously  recovered.  It  can  be  noted  that  the  A  matrix  is  11  x  11 
and  that  the  additional  equation  for  bin  60,10  was  not  included.  It  turns  out  that 
including  this  bin  does  not  add  any  additional  information. 
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B.  EXAMPLE 


Suppose  we  have  an  incoming  signal, 

s(t)  =  xo  +  x\  cos(27rt)  +  X2  cos(47rf)  4-  x%  cos(67rt)  +  X4  cos(87rf)  4- 
X5  COs(107T t)  +  Xq  COs(127 Tt)  +  Xj  COs(l47Tt)  +  X$  COS(167Tt)  + 
xq  cos(18^t)  +  2;  10  cos(207ri)  (66) 

All  xn  are  unknown.  The  signal  is  sampled  at  fs  1  =  10  Hz  and  /s2  =  11  Hz 
and  the  DFT  is  applied  to  both  with  iVi=10  and  A2=ll  points  respectively.  This 
produces  two  spectra  where  any  frequencies  from  6  to  10  Hz  are  aliased  into  bins  from 
0  to  5  Hz.  Both  DFT’s  are  normalized  and  the  results  for  bins  0  through  5  are  shown 
in  Figure  7.1a  for  /si=10  and  Figure  7.1b  for  fs 2=11. 

From  the  values  in  each  bin  it  is  possible  to  form  the  b  vector 

bT  =  [0.894  1.211  2.072  1.960  0.384  0.903  2.762  1.283  1.526  0.878  0.218]  (67) 

and  set  up  the  equation  x  =  A-1b  which  yields: 

xT  =  [0.218  0.047  0.679  0.751  1.441  0.384  0.519  1.321  0.532  0.847  0.831]  (68) 

The  solution  (68)  can  be  compared  with  the  actual  values  for  the  xn  which  is  done  in 
Figure  7.2. 

Figure  7.2a  shows  the  normalized  DFT  spectrum  that  results  from  sampling 
the  same  signal  at  20  Hz.  Figure  7.2b  plots  the  solution  (68).  The  solution  is  exact. 

The  solution  method  presented  so  far  does  little  good  if  only  signals  with  zero- 
phase  can  be  resolved.  In  the  next  chapter,  it  is  investigated  how  the  imaginary 
components  of  a  signal  alias  into  bins  and  the  results  are  used  to  generalize  the 
solution  for  signals  with  random  phase. 
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a) 


Bin 


b) 

Figure  7.1:  Normalized  DFT  for  frequencies  from  0  to  10  for  a)  f8 1  =  10 
and  b)  fs2  =  11. 
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Amplitude  Amplitude 


Frequency 


a) 


b) 

Figure  7.2:  Frequency  components  and  amplitudes:  a)  actual  and  b)  cal¬ 
culated 
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VIII.  RESOLVING  MULTIPLE  FREQUENCIES:  FIRST 

ALGORITHM 

A.  CONSIDERING  RANDOM  PHASE  SIGNALS 

Thus  far  attention  has  been  restricted  to  signals  with  zero  phase.  Conse¬ 
quently,  the  DFT  yields  an  all  real  spectrum.  When  the  frequency  components  have 
random  phase,  however,  a  modification  to  (65)  must  be  introduced. 

1.  Illustration  of  Random  Phase  by  Example 

Consider  a  signal  with  integer  frequency  components  from  0  to  10  Hz  given 
by: 

s(t)  =  x0  +  x\  cos(27 rt  4-  4> i)  +  x2  cos(4-7rt  4-  <t>2)  4 - b  zio  cos(207rf  +  0io)  (69) 


where  0n  represents  some  random  phase  for  each  frequency  component.  If  this  signal 
was  sampled  at  20  Hz,  DFT  applied  and  normalized,  bins  0  through  10  would  have 
the  values: 


Co 

Xq 

Cl 

x\  cos(0i)  4-  jxi  sin(<£i) 

C2 

= 

x2  cos(02)  +  jx 2  sin(02) 

(70) 

C10 

xio  cos(<f>io)  +  jx10  sin(0lo) 

where  c*,  represents  the  complex  value  of  bin  n  and  Ck  —  ak  +  jbk .  If  the  same  signal 
was  sampled  at  10  Hz,  the  DFT  applied  and  normalized,  the  real  part  of  each  bin 
will  be  given  by: 


al,10 

X\  COS(0i)  4-  Xg  COs(4>g) 

02,10 

x2  cos  (02)  +  x8  cos  (08) 

03,10 

= 

X3  COs(^3)  4-  x7  COs(^>7) 

O4,10 

X4  COS (^>4)  4-  Xq  cos(06) 

.  a5,10  . 

_  X5  cos(05) 
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2.  The  Imaginary  Component 

Consider  the  imaginary  component  of  each  bin.  If  the  DFT  is  treated  as  a  unit 
circle  as  in  Figure  3.1,  then  for  /s=10  Hz,  all  frequencies  below  5  Hz  lie  in  the  top 
half  of  the  unit  circle  and  appear,  without  aliasing,  in  bins  0  through  5  of  the  DFT. 
All  frequencies  from  6  to  10  Hz,  however,  lie  in  the  bottom  half  of  the  unit  circle  and 
alias  back  to  the  top  half  into  bins  0  through  5  of  the  DFT.  Since  the  real  part  of  the 
DFT  is  an  even  function,  the  real  components  of  the  signals  add  together  into  bins  0 
through  5.  The  imaginary  component  of  the  DFT  is  an  odd  function,  however,  and 
while  the  imaginary  component  of  frequencies  from  0  through  5  Hz  add  into  bins  0 
through  5,  the  imaginary  components  of  frequencies  from  6  through  10  Hz  subtract. 
In  other  words: 

frqio  x1sm(4>1)  -  xgsm(<t>9) 

62.10  *2  sin(<fo)  -  sin(^s) 

63.10  =  £3sin03)  -  x7sin (<£7)  (72) 

64.10  z4  sin (04)  -  xq  sin (<£6) 

65.10  J  L  *5  sin (<£5) 

B.  SOLVING  FOR  REAL  AND  IMAGINARY  COMPONENTS 

Based  on  how  the  imaginary  components  of  the  signals  alias  into  the  bins,  it 
is  again  possible  to  form  a  set  of  linear  equations. 


Z\  —  Zg 
Z2  —  Zg 
Z3  —  Z-j 

Z\  —  z§ 

z§ 

z§  —  z§ 

Z\  —  Z-j 
Z3  ~  Zg 
z2  —  Zg 

Z\  —  z10 

Zg 


(73) 
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■  Y0  1  rooooo  0  0 

Yi  11111—1—1 

r2  01111-1-1 

y3  ooiii-i-i 

r4  00011-1  o 

y5  =  00001  o  o 

Ye  0  0  0  0  1  -1  0 

Y7  0  0  0  1  1  -1  -1 

y8  ooiii-i-i 

Yg  0  1111-1-1 

_  Yw  J  1  1  1  1  1  1-1 


0  0  0  1  61,10 

—  1—1  0  0  62,10 

—  1  0  0  0  63,10 

0  0  0  0  64,10 

0  0  0  0  be, 10 

0  0  0  0  65,11  (75) 

0  0  0  0  64,11 

0  0  0  0  63,11 

—  1  0  0  0  62,11 

-1  -1  0  0  6i’n 

—  1  —1  —1  0  _  60,ii  . 
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which  solves  for  the  imaginary  component  while  the  real  component  is  solved  as  in 

(65): 


X0  ' 
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0 
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0 

“3,11 

X8 

0 
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-1 

-1 

-1 

1 

1 

1 

0 

0 

0 

“2.11 

Xg 

0 

-1 

-1 

-1 

-1 

1 

1 

1 

1 

0 

0 

ai.ii 

x10 . 

.  -1 

-1 

-1 

-1 

-1 

1 

1 

1 

1 

1 

0  . 

.  ao.n  . 

In  summary,  when  the  incoming  signal  has  components  with  random  phase,  both  the 
real  and  imaginary  components  may  be  solved  separately  by  the  equations  above. 


C.  THE  PIVOTAL  FREQUENCY 

An  inspection  of  the  specific  case  presented  illustrates  that  the  bin  information 
at  /  =  5  Hz  can  be  thought  of  as  pivotal  in  that  its  accuracy  affects  the  solution  for 
the  entire  spectrum.  Unfortunately,  while  the  phase  for  the  frequency  component  at 
/  =  5  Hz  (in  this  case)  is  random,  the  phase  of  Z$  is  not.  In  all  cases  for  an  N- point 
DFT,  where  N  is  even,  the  imaginary  component  of  the  Nj 2  bin  will  be  zero.  Thus, 
although  the  real  part  of  the  solution  will  be  accurate,  the  imaginary  component 
will  not.  In  order  to  solve  this  problem  with  the  DFT  for  the  simple  case  presented 
(fsi  =  10  Hz  and  fs 2  =  11  Hz),  the  signal  is  band-limited  from  0  to  5  Hz,  sample  at 
11  Hz,  apply  the  DFT  and  normalize.  The  value  of  this  DFT  at  bin  5  will  be  the 
correct  value  which  can  replace  the  values  in  (75)  and  (76)  for  65,10  and  a5)io. 


D.  EXAMPLE 

Suppose  there  is  an  incoming  signal, 

s(t )  =  0.218  +  0.047  cos(27rt  +  5. 1369)  +  0.679  cos(4tt*  +  4.7491)  + 
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(77) 


0.751  cos(6t ri  +  2.9044)  +  1.441  cos(87rt  +  5.9776)  + 

0.384  cos(107rt  +  3.9756)  +  0.519  cos(12tt*  4-  2.7604)  + 

1.321  cos(147rt  +  5.1817)  +  0.532  cos(16tt*  +  4.3290)  + 

0.847  cos(187rt  +  4.4121)  +  0.831  cos(20tt*  +  6.2024) 

Assuming  these  amplitudes  and  phases  are  unknown,  the  signal  is  sampled 
at  fsi  =  10  Hz  and  /s2  =  11  Hz,  the  DFT  is  applied  and  normalized.  Figure  8.1 
displays  the  real  and  imaginary  components  of  the  normalized  DFT  for  fs i  =  10  Hz. 
Figure  8.2  displays  the  real  and  imaginary  components  of  the  normalized  DFT  for 
fs2  =  11  Hz. 


0  12  3  4  5 


Bin 

b) 

Figure  8.1:  Normalized  DFT  components  at  fs i  =  10  Hz:  a)  real  compo¬ 
nent  and  b)  imaginary  component. 
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0  1  2  3  4  5 


Bin 

a) 


Bin 

b) 


Figure  8.2:  Normalized  DFT  components  at  fs 2  =  11  Hz:  a)  real  compo¬ 
nent  and  b)  imaginary  component. 

Note  from  Figure  8.1b  that  65,10*  =  0  which  is  the  imaginary  component  of 
the  fifth  bin  for  /sl=10.  This  is  not  the  actual  value  and  the  actual  value  must  be 
obtained  by  bandlimiting  the  signal  between  0  and  5  Hz,  sampling  again  at  11  Hz, 
and  applying/normalizing  the  DFT.  This  yields  the  correct  value  for  65,10  =  —0.2844. 
Now  we  may  form  the  two  vectors: 

aT  =  [-0.2312  -  0.1741  -  0.1325  -  0.8925  -  0.2580  -  0.7398  1.9717 

-0.9290  -  0.2256  0.8476  0.2180] 
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and 

bT  =  [0.7663  -  0.1852  1.3546  -  0.6266  -  0.2844  -  0.4775  0.7446 

0.6698  0.1306  0.0242  0.0000]  (78) 

Solving  X  =  A~la  and  Y  =  A^b  yields: 

'  Z0  ]  r  0.2180 
Zi  0.0194  -  j0.0428 

Z2  0.0249  -  j'0.6785 

Z3  —0.7300  +  ^0.1765 

Z4  1.3742  -  jf0.4335 

Z5  =  -0.2580  -  j 0.2844  (79) 

Z6  -0.4817  +  ^0.1931 

Z7  0.5975  -  j  1.1782 

—0.1990  —  j'0.4934 
Z9  -0.2505  -  jO.8091 

_  Zw  J  [  0.8283  -  j'0.0671  . 

A  simple  check  reveals  that  this  is  the  correct  result. 

This  simple  example  helps  to  illustrate  the  solution  method.  Systems  have 
been  simulated  in  which  /sl  =  1000  and  /s2  =  1001  and  exact  results  are  attained  in 
every  case.  One  of  the  drawbacks  of  the  solution  is  the  need  to  band-limit  the  signal 
and  sample  a  third  time  at  /s2  in  order  to  recover  the  phase  of  the  frequency  compo¬ 
nent  at  /si/2.  An  alternative  method  might  be  to  simply  not  transmit  that  frequency 
or  to  use  notch  filters  to  eliminate  that  frequency  and  the  problem  associated  with  it. 

E.  REDUCING  THE  COMPUTING  COST 

Concerning  the  computing  cost  of  the  algorithm,  if  two  full  matrix  multipli¬ 
cations  are  used  every  time  to  solve  for  the  spectrum,  computing  cost  becomes  a 
factor  as  each  matrix  multiplication  requires  0 (N2)  multiplications  and  additions. 
About  half  of  each  inverse  matrix  is  zeros  and  the  solution  could  be  hard-wired  more 
efficiently  to  reduce  flop  counts. 
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F.  GENERALIZING  THE  SOLUTION 


Having  presented  the  solution  method  by  example,  it  is  necessary  to  generalize 
the  result  for  any  two  sampling  frequencies. 


Theorem  8.1  Given  fs i  =  N\  where  Nj  is  an  even  number  and  fS2  =  N%  =  Ni  +  1, 
it  is  possible  to  form  two  N2  x  nonsingular  matrices  Ar  and  Ai  which  solve  the 


system  of  equations  ATX  =  a  and  AiY  = 
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b  given  by: 
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(81) 
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The  inverses  of  Ar  and  A*  exist  and  the  linear  system  may  be  solved  X  —  Arla 


and  Y  =  Ai1b  given  by: 
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(83) 


All  methods  for  solving  undersampled  signals  so  far  have  required  at  least  two 
sampling  frequencies.  The  additional  complexity  is  a  trade-off  for  the  additional  infor¬ 
mation  gained.  In  the  next  chapter,  however,  a  second  algorithm  is  introduced  which 
requires  only  one  sampling  frequency  and  simplifies  the  overall  solution  considerably. 
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IX.  RESOLVING  MULTIPLE  FREQUENCIES:  SECOND 

ALGORITHM 


A.  BACKGROUND 

So  far  the  solution  method  has  required  that  an  undersampled  signal  be  sam¬ 
pled  more  than  once  with  a  different  sampling  frequency  each  time.  It  is,  however, 
possible  to  undersample  a  signal  more  than  once  with  a  single  sampling  frequency 
yet  still  recover  the  original  spectrum.  In  fact,  the  groundwork  for  such  a  system  has 
already  been  laid  in  Chapter  VIII.  Recalling  the  discussion  in  Chapter  VIII,  Section 
C,  the  problems  encountered  with  the  pivotal  frequency  required  that  the  signal  be 
band-limited  between  0  and  41  to  recover  the  actual  value  of  bin  c/^  ,  .  It  appeared 
to  be  a  waste  of  time  but  a  necessary  evil  to  recover  the  correct  bin  information  for 
this  pivotal  bin.  Options  discussed  to  remove  this  apparent  constraint  were  to  not 
transmit  the  frequency  or  use  a  notch  filter  to  eliminate  it  before  processing. 

B.  USING  THE  BANDLIMITED  SPECTRUM 

Consider  a  signal  s(t)  which  contains  frequency  components  from  0  to  /  Hz 
of  unknown  magnitude  and  phase.  Previously,  the  solution  would  have  required  that 
s(t)  be  sampled  at  both  fs i  =  /  Hz  and  /s2  =  /  +  1Hz  in  order  to  recover  the  original 
spectrum.  As  discussed,  the  solution  would  also  require  that  s(t )  be  additionally 
bandlimited  between  0  and  // 2  Hz  in  order  to  recover  the  correct  value  for  cQjsl- 
Discarding  all  information  gained  from  the  band-limited  spectrum,  (82)  and  (83) 
would  then  be  used  to  solve  for  the  original  spectrum.  In  fact,  the  solution  method 
wastes  time  because  the  band-limited  spectrum  has  no  aliased  components  and  the 
bin  values  correspond  to  the  actual  values  for  Zq,  Z\, ...,  Z i . 
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1. 


Recovering  the  First  Half  of  the  Spectrum 


Consider  a  signal,  s(t)  containing  frequency  components  from  0  to  10  Hz.  If 
the  signal  is  first  bandlimited  between  0  and  5  Hz  and  then  sampled  at  11  Hz,  the 


actual  frequency  coefficients  are  given  by: 


'  ■ 

C0,11  L 

Zi 

Cl,ll  L 

z2 

c2,llt 

z3 

c3,lli 

z4 

c4,11l 

z5 . 

.  c5,11l  . 

(84) 


where  Zn  =  Xn  +  jYn  represents  the  actual  real  and  imaginary  component  at  the 
frequency  of  n  Hz  and  cnjsL  =  anjsL  +  jbnjsL  corresponds  to  the  value  of  the  nth  bin 
at  a  sampling  frequency  of  fs.  The  additional  L  subscript  denotes  that  this  is  the 
band-limited  version  of  the  signal. 


2.  Recovering  the  Other  Half  of  the  Spectrum 


Suppose  that  on  another  channel  the  signal  s(t )  is  not  band-limited.  If  s(t)  is 
sampled  at  11  Hz  and  the  DFT  is  applied  and  normalized  the  real  components  will 


go  into  bins  as: 


a0,ll 
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<*1,11 

X\  4-  Xio 

<*2,11 

x2  +  x9 

<*3,11 

x2  +  x$ 

<*4,11 

x4  +  x7 

.  <*5,11  . 

x5  +  x6  _ 

while  the  imaginary  components  go  into  bins  as: 


bo,n 

01 

bi,n 

o 

i 

b2,u 

T2-y9 

&3,11 

Y3-Ys 

^4,11 

y4-y7 

.  h,n  . 

Y5-Y6 

(85) 


(86) 
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Substituting  the  values  for  the  Zn' s  from  (84),  (85)  and  (86)  may  be  re-written 


as: 


'X10 ' 

<*1,11  —  *1,1U 

*9 

a2,ll  ~  a2,llL 

= 

*3,11  —  *3, llz, 

a7 

*4,11  “  *4, llL 

a6 

_  *5,11  ~  *5, 11L 

'  Yio  ‘ 

51,11  —  &1,11L 

fa, 11  -  fa, 111 

F8 

=  — 

53,11  -  fa, 11l 

r? 

54,11  —  fa, 11 L 

n 

_  fa, 11  -  fa,llL  . 

Combining  (84),  (87),  and  (88)  yields: 


'  Zo  - 

*0,lli, 

Zi 

*l,llz 

z2 

*2,11l 

z3 

*3,11l 

Z4 

*4,11  L 

Z5 

= 

*5, 11L 

Zq 

(  *5,11  —  *5,11!,  ^ 

Z7 

*4,11  —  *4,11  L 

Zg 

conj 

*3,11  “  *3, Hi, 

Z9 

*2,11  -  *2,11!, 

.  Z10  . 

. 

*1,11  ~  *1,11l  / 

(87) 


(88) 


(89) 


As  can  be  seen,  the  solution  is  much  simpler  requiring  only  one  sampling  frequency 


and  far  fewer  computations. 


C.  GENERAL  SOLUTION 

Figure  9.1  displays  a  possible  architecture  to  resolve  an  undersampled  signal 
using  the  method  described.  The  incoming  signal  is  split  into  two  branches.  The 
top  branch  is  the  band-limited  branch  and  solves  for  the  first  half  of  the  spectrum. 
This  information  is  then  combined  with  the  information  from  the  bottom  branch  to 
recover  the  undersampled  spectrum. 

The  solution  (89)  can  be  generalized: 
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Theorem  9.1  Let  the  sampling  frequency  of  a  system  be  fs  where  fs  is  any  odd  num¬ 
ber.  A  spectrum  containing  frequencies  between  0  and  fs  —  1  Hz  may  be  unambiguously 


recovered  provided  that  a  copy  of  the  signal  is  also  band-limited  to  recover  the  spectrum 


between  0  and  Ap-. 


This  information  is  then  used  to  recover  the  other  half  of  the 


spectrum  as  in: 


Zo 

°0  JsL 
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Zi 

CUsL 

z2 

C%fsL 

clsf±-2JsL 

Z^-l 

C'kfi-1  ,fsL 

Z  fs-l 

2 

'  ~,fs  °^JsL  ' 

C  ^2  1  ~li  fs  C  'ky-  —IJsL 

Z^+3 

chfl-2,fs  ~  C^i-2 ,fsL 

conj 

Z  fs— 3 

c3,/s  -  c3  JsL 

Z  fs— 2 

C2Js  -  C2 JsL 

Zfs- 1  J 

\  Citfs  -  cijsL  y 

where  Zn  —  Xn  +  j  Yn  and  cn  m  —  an  m  +  jbn  rn . 


(90) 


It  would  be  nice  to  know  that  (9.1)  puts  the  Nyquist  criterion  to  rest  once  and 
for  all  but  a  fair  treatment  requires  that  the  problems  with  the  methods  described 
be  explored.  The  next  chapter  discusses  the  current  limitations  of  the  algorithms  as 
well  as  potentially  resolving  the  same  limitations. 
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Figure  9.1:  Channel  architecture  to  resolve  an  undersampled  signal  utiliz¬ 
ing  a  single  sampling  frequency 
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X.  PROBLEMS  WITH  THE  SOLUTION 


A.  SIGNALS  THAT  CAUSE  DFT  LEAKAGE 

So  far  only  frequencies  which  directly  resolve  into  single  bins  without  leakage 
have  been  considered.  There  is  a  problem  with  the  solution  when  the  incoming 
frequency  components  lie  in  between  DFT  bins  and  leakage  occurs.  When  the  DFT 
is  applied  to  a  sequence  with  frequencies  that  correspond  directly  to  DFT  bins,  the 
result  is  that  the  components  resolve  all  of  their  energy  into  the  bins  of  the  DFT. 
This  is  not  so  with  frequencies  that  “straddle”  bins. 

B.  DFT  LEAKAGE  EXAMPLES 

Consider  a  signal  containing  just  one  frequency  at  2.3153  Hz,  given  by: 

s(t)  =  cos(27r(2.3153)t)  (91) 

Note  that  the  signal  has  zero  phase.  Unfortunately,  instead  of  resolving  into  just  bin 
2  when  sampled,  the  side-lobes  of  the  DFT  distribute  a  portion  of  the  energy  into 
many  bins. 

1.  Algorithm  of  Theorem  8.1 

Using  the  algorithm  of  Theorem  8.1,  the  signal  is  sampled  at  both  10  and  11 
Hz,  the  DFT  and  normalized  with  results  shown  in  Figure  10.1.  Displayed  are  the 
real  and  imaginary  components  of  the  DFT.  These  DFT  results  along  with  (82)  and 
(83)  yield  the  spectrum  shown  in  Figure  10.2a  for  the  real  part  and  Figure  10.2b  for 
the  imaginary  part.  The  result  is  clearly  inaccurate. 

2.  Algorithm  of  Theorem  9.1 

The  algorithm  of  Theorem  9.1  is  not  nearly  as  inaccurate  for  a  single  frequency. 
Figure  10.1c  and  Figure  lO.ld  correspond  to  both  the  band-limited  and  normal  spec- 
trums  which  would  be  used  to  solve  for  the  original  spectrum.  The  answer  would  be 
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exact.  The  algorithm’s  strength  lies  in  the  fact  that  half  of  the  solution  is  determined 
exactly  in  the  top  branch  of  Figure  9.1.  The  solution  is  not  totally  immune  to  the 
effects  of  DFT  leakage  however.  As  an  example,  assume  a  signal  s(t)  has  frequency 
components  in  \  Hz  increments  from  0  to  10  Hz  (0,  .5,  1,  1.5,  ...,  9.5,  10).  The 
sampling  frequency  is  11  Hz.  Figures  10.3a  and  10.3b  display  the  real  and  imaginary 
spectra  of  the  band-limited  version  of  the  signal  after  sampling  and  DFT  normaliza¬ 
tion  while  figures  10.3c  and  10.3d  display  the  non-bandlimited  spectra  after  sampling 
and  DFT  normalization.  Utilizing  (90),  the  original  spectrum  is  calculated  and  plot¬ 
ted  in  Figures  10.4a  and  10.4b  while  Figures  10.4c  and  10. 4d  show  the  actual  values. 
The  solution  is  inaccurate. 

C.  LEAKAGE  AND  SOLUTION  INSTABILITY 

It  is  not  too  difficult  to  discern  why  the  solution  methods  are  unstable.  One 
might  think  that  perhaps  the  leakage  from  all  the  other  frequencies  would  somehow 
combine  in  such  a  way  that  the  relative  magnitudes  of  the  resulting  spectrum  would 
remain  intact.  Unfortunately,  the  symmetrical  number  system  and  the  linear  solutions 
presented  both  depend  on  nearly  all  energy  being  within  the  bin  that  it  is  expected 
to  be  in.  As  an  example  consider  if  an  incoming  frequency  at  2  Hz  correctly  processes 
totally  into  bin  2  at  10  Hz  but  somehow  processes  into  bin  3  at  11  Hz.  Whereas  the 
correct  solution  is  for  all  energy  to  be  at  f  =  2  Hz  in  the  final  spectrum,  this  small 
bin  error  will  result  in  calculating  that  all  energy  belongs  to  f  =  8  Hz. 

D.  SOLVING  THE  PROBLEM? 

What  is  needed  ideally  is  some  sort  of  transformation  that  places  all  energy 
in  the  expected  bins  without  leakage.  For  instance,  for  fs\  =  10  it  is  desired  that 
/  >  9.5  and  0  <  /  <  0.5  to  resolve  into  bin  0,  0.5  <  /  <  1.5  to  resolve  into  bin  1, 
•  ■  and  8.5  <  /  <  9.5  to  resolve  totally  into  bin  9.  Clearly  no  such  transformation 
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exists  which  will  resolve  a  spectrum  exactly  without  leakage.  There  is  one  direction, 
however,  that  may  contain  some  promise. 

Consider  again  our  signal,  s(t),  with  a  single  frequency  at  2.3153  Hz.  Sample 
again  at  10  and  11  Hz  but  this  time  sample  for  100  seconds  instead  of  just  1  second. 
Applying  the  DFT  and  normalizing,  the  DFT  now  has  a  bin  resolution  of  0.01  seconds 
instead  of  just  1  second.  This  means  that  any  incoming  frequency  at  any  integer 
multiple  of  0.01  Hz  will  resolve  exactly  into  a  single  bin  while  all  others  will  cause 
leakage.  The  results  of  the  DFTs  are  shown  in  Figure  10.5.  Figure  10.5  shows  the 
effect  of  increasing  the  bin  resolution  which  has  the  effect  of  containing  the  energy 
within  the  area  of  interest.  The  sidelobes  decay  well  before  they  “leak”  much  into  any 
of  the  other  bin  boundaries.  Still,  while  the  energy  is  now  contained  in  the  region  of 
interest,  there  needs  to  be  a  method  to  recover  the  overall  real  and  complex  portions 
of  the  signal  that  belong  to  each  integer  bin.  In  the  case  studied,  the  goal  is  to  find  the 
energy  that  resolves  into  the  five  bins  for  fs i  =  10  Hz  and  fs 2  =  11  Hz,  respectively, 
and  use  the  linear  equations  presented  to  solve  for  the  spectrum  with  a  bin  resolution 
of  1  Hz. 

For  the  algorithm  of  Theorem  9.1  the  same  benefits  apply.  In  fact,  increasing 
the  bin  resolution  may  hold  even  more  promise.  The  algorithm  holds  up  fairly  well 
when  only  a  little  leakage  occurs  but  falls  apart  when  leakage  increases. 

It  should  be  emphasized  that  10  and  11  Hz  have  been  used  as  sampling  fre¬ 
quencies  throughout  for  simplicity  of  presentation.  Preliminary  study  has  shown 
that  a  separation  of  1  between  the  two  sampling  frequencies  gives  the  best  dynamic 
range  with  the  lowest  sampling  frequencies.  It  can  be  shown  easily,  for  instance,  that 
the  dynamic  range  of  the  algorithm  of  theorem  8.1  is  still  0  to  10  Hz  for  sampling 
frequencies  of  10  and  12  Hz. 
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Another  aspect  of  theorem  8.1  is  that,  unlike  the  multiplicative  nature  of 
the  dynamic  range  with  the  symmetrical  number  system,  the  dynamic  range  of  the 
solution  is  additive.  It  depends  upon  the  number  of  linearly  independent  equations 
that  can  be  set  up  based  upon  the  number  of  bins  for  each  sampling  frequency  used. 
For  this  reason,  this  algorithm  would  be  most  cost-performance  efficient  using  just 
two  sampling  frequencies  as  presented. 
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Figure  10.1:  Normalized  DFT  of  sinusoid  at  /  =  2.3153  Hz:  a)  real  com¬ 
ponent  for  fsi  =  10  Hz,  b)  imaginary  component  for  fs i  =  10  Hz,  c)  real 
component  for  fs 2  =  11  Hz,  and  d)  imaginary  component  for  fs 2  =  11  Hz. 


a) 


b) 

Figure  10.2:  Spectrum  calculated  for  single  frequency  at  2.3153  Hz:  a) 
real  component  and  b)  imaginary  component 
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Figure  10.3:  Spectra  of  a  signal  with  frequency  components  from  0  to 
10  Hz  in  1/2  Hz  increments:  a)  real  component  of  sampled  bandlimited 
signal,  b)  imaginary  component  of  bandlimited  signal,  c)  real  component 
of  sampled  original  signal,  and  d)  imaginary  component  of  original  signal 
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Frequency  (Hz) 

4) 

Figure  10.4:  Spectra  of  a  signal  with  frequency  components  from  0  to 
10  Hz  in  1/2  Hz  increments:  a)  calculated  real  component,  b)  calculated 
imaginary  component,  c)  actual  real  component,  and  d)  actual  imaginary 
component. 
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Figure  10.5:  Normalized  DFT  with  0.01  Hz  resolution  of  sinusoid  at 
/  =  2.3153  Hz:  a)  real  component  fo  fsl  =  10  Hz,  b)  imaginary  compo¬ 
nent  for  fsi  =  10  Hz,  c)  real  component  for  fs 2  =  11  Hz,  and  d)  imaginary 
component  for  fs 2  =  11  Hz. 
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XI.  CONCLUDING  REMARKS 


The  essential  contribution  of  this  thesis  is  detailing  algorithms  or  solution 
methods  that  may  hold  the  key  to  recovering  an  undersampled  spectrum  using  one  or 
more  sampling  frequencies.  Using  the  Symmetrical  Number  System  and  the  Discrete 
Fourier  Transform  it  is  a  very  simple  process  to  resolve  one  or  two  undersampled 
frequencies  within  the  dynamic  range  of  the  SNS.  Additionally,  when  only  resolving 
one  or  two  frequencies,  each  pairwise  relatively  prime  sampling  frequency  multiplica- 
tively  extends  the  dynamic  range  of  the  system.  When  resolving  an  entire  spectrum 
of  frequencies  the  solution  becomes  more  complicated  but  it  is  still  a  straightforward 
process.  The  number  of  unique  bins  for  each  sampling  frequency  dictates  the  number 
of  linearly  independent  equations  that  can  be  formed  and  hence  dictates  the  dynamic 
range.  Using  at  most  two  sampling  frequencies  demonstrated  it  has  been  shown  that 
the  spectral  width  that  can  be  recovered  using  either  method  is  approximately  double 
the  width  predicted  by  Nyquist. 

Future  efforts  will  attempt  to  resolve  two  important  questions.  First,  is  it 
possible  to  apply  the  solution  to  a  signal  which  has  frequency  components  that  don’t 
conveniently  lie  directly  at  bin  frequencies?  The  answer  to  this  question  may  lie  in 
increasing  bin  resolution  and  then  using  a  mathematical  method  to  determine  fairly 
accurately  the  total  energy  that  actually  belongs  to  each  integer  frequency.  The 
second  question:  What  happens  when  more  than  two  sampling  frequencies  are  used? 
The  matrix  solution  for  two  sampling  frequencies  is  very  simple  and  the  solution 
matrix  has  nice  properties.  That  may  or  may  not  be  the  case  when  three  or  more 
sampling  frequencies  are  used. 

In  spite  of  problems  that  remain  with  the  methods  presented,  the  thesis  is 
significant  as  it  presents  a  direction  for  resolving  an  undersampled  spectrum.  Provided 
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that  the  leakage  problem  is  solved,  either  or  both  of  the  methods  may  hold  the  key 
to  greatly  extending  the  ability  to  sample  and  resolve  extremely  high  frequencies  at 
relatively  low  sampling  rates. 
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